Critical exponents of the anisotropic Bak-Sneppen model 
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We analyze the behavior of spatially anisotropic Bak-Sneppen model. We demonstrate that a 
nontrivial relation between critical exponents r and /i = d/D, recently derived for the isotropic 
Bak-Sneppen model, holds for its anisotropic version as well. For one-dimensional anisotropic Bak- 
Sneppen model we derive a novel exact equation for the distribution of avalanche spatial sizes, and 
extract the value 7 = 2 for one of the critical exponents of the model. Other critical exponents are 
then determined from previously known exponent relations. Our results are in excellent agreement 
with Monte Carlo simulations of the model as well as with direct numerical integration of the new 
equation. 



Ever since its introduction five years ago, the Bak- 
bneppen (BS) model § was a subject of considerable 
theoretical interest. Its relevance relies on the fact that 
it provides a very simple mechanism of Self-Organized 
Criticality 0. In fact the Bak-Sneppen model is the sim- 
plest representative of a broad class of extremal models, 
which all naturally evolve towards a scale-free stationary 
state ||. An extremely rich dynamic critical behavior 
arising out of truly minimalistic dynamical rules has in- 
spired numerous analytical and numerical investigations 
of the Bak-Sneppen model [p|-^0[. 

In this work we introduce and study an anisotropic 
version of the original Bak-Sneppen model. In one di- 
mension the dynamics is as follows: the configuration of 
the system is fully defined by the value of variable fi for 
each lattice site i. At every time step the smallest vari- 
able in the system and that at its right nearest neighbor 
are replaced with new random numbers independently 
drawn from the distribution V{f) = p[ . Contrary 
to the isotropic BS model, where both nearest neighbor 
variables are updated, only the variable in the preferred 
direction is updated here. Other mechanisms of intro- 
ducing spatial anisotropy to the rules of the original BS 
model were recently studied in Refs. 12 13||. As we shall 
see, the universality of the critical behavior manifests it- 
self in the fact that any realization of the anisotropy in 
the original BS model gives rise to the same set of crit- 
ical exponents Jl|] . The generalization of our version of 
anisotropic BS model to higher dimensions is straight- 
forward: only d neighbors of the global minimal site lo- 
cated in positive directions of corresponding coordinates 
are updated along with it. 

This work is devoted to analytical and numerical study 
of exponents of the anisotropic Bak-Sneppen model. The 
main observations used in the analytical part of this 
study are: i) the general scaling theory of Ref. [|| devel- 
oped for an arbitrary extremal model reduces the num- 
ber of independent critical exponents to just two; ii) 
the relation t(jx) between two remaining exponents r 



and \x = d/D, recently derived in |8j,|9| for the isotropic 
Bak-Sneppen model in an arbitrary dimension, holds for 
its anisotropic version as well. Finally, Hi) in the one- 
dimensional anisotropic BS model the exponent 7 = 2, 
describing the divergence of the average avalanche size, 
can be derived from the exact master equation for the 
probability distribution of avalanche spatial sizes. This 
equation, which we derive in this paper, compliments 
that for the probability distribution of avalanche tem- 
poral durations From the exact result 7 = 2 it fol- 
lows that the values of exponents r and [i in the one- 
dimensional anisotropic BS model are related via r = 2/i. 
Using the approximate form of the function r(/x), de- 
rived in ||, we can give an analytical estimate of the 
exponents r and \x in ID anisotropic BS model. Indeed, 
they must lie at the intersection point of the r = 2/i 
line and the r(/x) curve, valid for an arbitrary BS model 
H . These analytical estimates are in excellent agreement 
with the direct numerical integration of the exact master 
equation for PDF of the avalanche spatial sizes, giving 
H = 0.588(1) and t = 1.176(2), as well as Monte Carlo 
simulations performed by us and by he authors of Refs. 

US 

The scale-free stationary state of an arbitrary extremal 
model can be characterized by a number of critical expo- 
nents. Several scaling relations, reviewed and discussed 
in Ref. reduce this variety to only two independent 
exponents such as r for the power law in the distribu- 
tion P(s) ~ s~ T of avalanche temporal durations s, and 
the dynamic exponent fi. The latter exponent relates 
avalanche temporal duration s to its spatial volume V(s) 
as V(s) ~ s' 1 . Here the spatial volume is defined as 
the number of distinct sites updated at least once during 
this avalanche. In the notation of Ref. ||] /j = d/D (or 
/i = d cov /D if the set of updated sites is not compact, 
but instead forms a fractal of dimension d cov ). 

The Bak-Sneppen model in an arbitrary dimension and 
with an arbitrary anisotropy has an additional simplifica- 
tion 0] : since variables are simply replaced with new ran- 
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dom numbers and have no memory about their previous 
values, the dynamics within a single avalanche is totally 
independent from what happened before it started. This 
observation j| enables one to simulate /-avalanches (for 
definitions see ||) for an arbitrary value of /, which can 
be above as well as below the critical point, without spec- 
ifying variables at passive sites (those with /j > /) prior 
to avalanche. Another important consequence of the ab- 
sence of inter-avalanche memory is that in any variant 
of BS model there are no correlations between sizes of 
subsequent avalanches. This statement is not a result 
of any kind of mean-field approximation. Rather it is a 
clear logical consequence of the dynamical rules of the 
Bak-Sneppen model. 

In @] one of us derived a master equation for the distri- 
bution P(s, /) of avalanche durations s, valid for a gen- 
eral BS model. Let us recall briefly the sequence of logical 
steps leading to this equation. The starting point is the 
analysis of the signal (minimal number as a function of 
time) /min(i) of the model using an auxiliary parame- 
ter /. The intersection of this signal with the horizontal 
line drawn at / identifies the sequence of /-avalanches, 
following one another. In other words, if f m i n (t) > f, 
fmin(t + k) < f for 1 < k < s, and f min (t + s) > f, 
we say that an /-avalanche of size (temporal duration) s 
has occurred |9j. The sequence of /-avalanches is char- 
acterized by the distribution P(s, /) of their sizes s. One 
can investigate how this distribution changes under an 
infinitesimal increase of / from / to f + df. The change 
occurs simply because when the horizontal line at / is 
lifted, some intersections with the signal f m i n (t) disap- 
pear. This means that two consecutive /-avalanches of 
(temporal) sizes Si and S2 merge into a single / + df- 
avalanche of size s\ + S2- The occurrence of this event im- 
plies that at least one of the V\ ~ sites, updated during 
the course of the first avalanche, has f < fi < f + df. 
Taking into account that as we argued above subsequent 
avalanches in Bak-Sneppen model arc uncorrelated, one 
can write the balance of loss and gain of /-avalanches of 
size s as / is increased. The resulting master equation 
for the distribution P(s, f) is given by 

s-l 

dfP(8, /) = -S*P(s, /)+£ S»P( Sl , f)P(s - 8 U /). 

si=l 

(1) 

Strictly speaking in order for the above equation to be 
exact one needs to replace s M with the average number 
of updated sites V(s, /), where the average is taken over 
all /-avalanches of size s. As was observed numerically 
this quantity has an insignificant /-dependence and its 
large s asymptotics is well described by the power law 
V(s,f) ~ As 1 - 1 . It was suggested in || and later con- 
vincingly confirmed numerically in |5) that the critical 
exponent \x uniquely determines the scaling properties of 



P(s, /). This justifies our substitution of V(s,f) by its 
asymptotical form s^ in Eq. ([[]) (the constant A in front 
of s M was absorbed by redefinition of /). 

The numerical integration of equation (|l|) with initial 
conditions P(s, f = 0) = S S: i shows that as / approaches 
some critical value f c (n), P{s,f) develops a power law 
form with a diverging cutoff: P(s,f) — s~ T F(s a (/ c — /)). 
Above f c there is a finite probability Poo(f) to start 
an avalanche that never ends. The possibility of such 
event shows up in Eq. ([!]) through the "normalization 
catastrophe", when -P( s ; /) f° r f > fc starts to 

fall below unity. This deviation is attributed to the 
appearance of the "infinite avalanche" with probability 
Poo(f) ~ (/ — fc) 13 ■ This way the overall normalization 
J2^Li -P( s i /) +Poo(/) = 1 is satisfied at all /. The prop- 
erties of Eq. (Q) depend on the critical exponent /i. Given 
the value of this dynamic critical exponent, the remain- 
ing exponents r, a = [i + 1 — r, 7 = (2 — r)/a, and 
/3 = (t — l)/cr, as well as the scaling function F{x) of 
the avalanche distribution near the critical point follow 
from Eq. (|l|). In |l the function r(/i) and the scaling 
form F{x) were determined by numerical integration of 
Eq. (Q) and the expansion around the mean-field point 
/i = 1, t = 3/2. It was shown H that to a second order 
in 1 — fx the function t(/i) is given by 

t(h) = 1-5 - (1 - ») + c(l - tif + 0((1 - A 4 ) 3 ); 
c=-(7 e +ln2-l)~ 0.3605. (2) 

Here j e ~ 0.5772 is the Euler's constant. 

Our new results for the one-dimensional anisotropic 
BS model are based on the exact equation for the prob- 
ability distribution function Q(r, /) of spatial sizes r of 
/-avalanches. The derivation of this equation is very sim- 
ilar to the derivation of (]l|), briefly outlined above. The 
main difference lies in the fact that when an avalanche 
of spatial size n merges with that of size ri they can 
overlap to form any spatial size between max(ri, r^) and 
7*1 +T2 — 1. Contrary to this, in the temporal domain the 
merging of avalanches of temporal durations Si and S2 al- 
ways produces an avalanche of temporal duration Sj +S2- 
Let us analyze how Q(r, /) changes when / is increased 
by an infinitesimal amount df. Some avalanches of size r 
will merge with the next avalanche. This event can only 
occur if one of the r sites, updated in the first avalanche, 
happens to host the smallest number right after the 
avalanche is finished. Each of these r sites has a variable 
fi, which was randomly drawn from T'(f) — during 
the course of the avalanche. At the end of this avalanche 
by the very definition of an /-avalanche all these r sites 
have fi > f. We can therefore regard the /j's on these 
sites as randomly drawn from an exponential distribu- 
tion normalized between / and 00. The probability that 
a particular fi is in the interval [/, / + df] , to linear order 
in df, is just df. The probability that at least one of the 
r sites has an fi in the interval [/, / + df] is rdf + 0(df 2 ). 
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This implies that the number of /-avalanches which will 
merge with the next one when / is raised to / + df is 
dQ(r)\i oss = —rdfQ(r, f) + 0(df 2 ). Let us now consider 
a merging event between two /-avalanches of size n and 
r 2 resulting in an /+d/-avalanche of size r. There are two 
scenarios of how this can happen: (i) the rightmost point 
of the second avalanche is at distance r from the leftmost 
(starting) point of the first avalanche. The constraint on 
possible values of r± and r 2 imposed by this scenario is 
max(ri, r 2 ) < r < r\ + r 2 — 1; (ii) r\ = r, and the second 
avalanche is fully contained within the first one. In the 
former case, the values of r, r±, and r 2 uniquely specify 
the initial site of the second avalanche. Therefore, the 
probability of this event to occur is just the probability 
df + 0(df 2 ) that this site has /, G [/, f + df]. However, in 
the latter case the starting point of the second avalanche 
can be any of the first r — r 2 sites of the first avalanche. 
This event occurs with a probability (r — r 2 )df + 0(df 2 ). 
Putting all these terms together we find: 



d f Q(r,f) = -rQ(r,f)+ j^Q(n,/) £ 

r 

+ Q(rJ) ^2(r-r 2 )Q(r 2 ,f). (3) 

r 2 = l 

This is an exact equation for the distribution of spa- 
tial sizes of /-avalanches. Unlike our previous results, 
its validity does not require any scaling assumptions, 
such as V(s, /) ~ s M used in the derivation of Eq. 
(|l|). The equation (||) has to be solved with the ini- 
tial condition Q(r, f = 0) = <5 r ,2- Indeed, in the one- 
dimensional anisotropic BS model at any time step (or 
/ = 0-avalanchc for that matter) r = 2 sites are updated. 

A similar but more complicated equation can be 
written for the isotropic one-dimensional Bak-Sneppen 
model. The basic object of this equation is the proba- 
bility distribution Q(ri,r 2 , /) of /-avalanches replacing 
precisely r\ sites to the left of the starting point, and r 2 
sites to the right of it. The initial condition is given by 
Q( r i7 r 2,0) = cVi,i<5r 2 ,i- Although we were unable to de- 
rive any analytical exponent relations from this equation, 
it should be stressed that it contains all properties of the 
one-dimensional isotropic BS model, and in principle its 
numerical integration constitutes a viable alternative to 
Monte Carlo simulations of the model. 

Yet another variant of Eq. (||) can be written for 
the anisotropic Bak-Sneppen model in dimensions higher 
than one. Let r to denote the spatial extent of the 
avalanche, measured along the diagonal (1, 1, . . . , 1) of 
the d-dimensional space. In projection to this axis the 
starting point of an avalanche is always the leftmost point 
in the avalanche. In order to describe the shape of the 
region covered by avalanche one needs to introduce the 
new exponent defined by the average number of up- 
dated sites n pro j(r) ~ r*> projected onto the same point 



at the distance r along the diagonal from the starting 
point. Then the total number of points covered by an 
avalanche of size r, n cov (r) 

anisotropic BS model n pro j(r) = 1, and, therefore, £ = 0. 
By retracing arguments that led to the derivation of Eq. 
(0) it is easy to see that its higher-dimensional version 
can be written as 

d f Q(rJ) = -r 1+ ^Q(rJ) 

r r 

+ $>iQ( r i>/) £ Q(r a ,/) 

ri — 1 T2— r— ri+1 

r r—r s — l 

+ Q(r,f)Y,r C s £ Q(r 2 ,f). (4) 

r s — 1 T2 — 1 

The drawback of this equation is that, similar to Eq. (Q), 
it requires the input of an additional parameter (critical 
exponent) £. 

Note that the Equation (^) for Q(r, f) involves only 
Q(r',f) for r' < r. Therefore, in principle this distri- 
bution can be computed numerically for r < R to the 
desired accuracy. Forward numerical integration of Eq. 
(||) shows that as / — > f c ~ 1.2865, Q(r, /) develops 
a power law behavior with an exponent r r = 1.299(3) 
(see Fig. 1). It is easy to see that in one dimen- 
sion the exponent r r of the distribution of avalanche 
spatial sizes is related to the more familiar exponent r 
of the distribution of their temporal durations through 
r r = (t — l)//x+l. Indeed, since asymptotically r — As 11 , 



P(s > s ) = P(r > As%) ~ (As%) 
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ing powers of so in this expression one gets the above 
exponent relation. To find values of critical exponents 
hidden in Eq. (|^) we study the behavior of moments 
(r n ) of the distribution Q(r, /) as a function of /. These 
satisfy the equation: 

oo r + 



df(r n ) = -{r n+1 )+J2 E 



7*+ = l 7*_=1 



ro. +t*_ —1 



2 E P n + rl(r + -r_) 



p=r + 



Q(r +) /)Q(r_,/). (5) 



where the sum over r\ and r 2 has been trans- 
formed into a sum over r + = max(ri,r2) and 
r_ = rain (ri, 7*2). For n = 1 the Eq. (|^) reads 



df{r) 



(r 2 ) + iE;:=i( 



r_ + r_ra 



r_)Q(r+,/)Q(r_,/) - -<r 2 ) + E^i E~=i VI + A + 
rir 2 - min(ri,r 2 )]/2 Q(n, f)Q(r 2 , /). For / < f c , when 
there are no infinite avalanches and the avalanche distri- 
bution Q(r, /) is normalized to unity, one gets 



df(r) 



(min(ri, r 2 )) 



(6) 



Close to the critical point we can neglect the term 
(min(ri, r 2 )), since it diverges slower than (r) 2 and we 
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are left with df(r) 
tion 

(r) 



(r) 2 /2. This has an obvious solu- 



fc-f 



O 



1 



(fc-ff 



(7) 



It has been known for some time that in the Bak-Sneppen 
model as well as in several other extremal models the 
amplitude of the divergence of (r) as / — > /~ is given 
by the critical exponent 7. This exponent describes the 
critical behavior of the average avalanche duration as 
(s) ~ |/c — /|~ 7 . For the original derivation of this fact 
[3 see Eq. (17) in ||. Actually, as it was shown in 
[^7 this fact can be also derived from Eq. (|lj). Indeed, 
from this equation it is easy to see that in the absence 
of infinite avalanche the first moment of P(s, f) obeys 
df(s) = (s fl )(s). Therefore, in the critical region one has 
(r) = = d f (ln(s)) = 7 /(/ c - /) + 0((/ c - /)- 2 )! 
From Eq. (Q) we conclude that in the anisotropic one- 
dimensional Bak-Sneppen model 



7 = 2. 



(8) 



Using the scaling relation || 7 = (2 — r)/a = (2 — 
t)/(1 + /j, — r), we readily find that Eq. (||) implies 
t = 2/i. This, combined with the r(/x) relation found 
in ||, gives jU = 0.58(1) and r = 1.16(2) as coordinates 
of the intersection point (see Fig. 2). The uncertainty in 
these numbers comes from our approximate knowledge 
of the systematic errors in the t(/z) curve, determined by 
numerical integration of Eq. (Q). Presently, this integra- 
tion was performed and approximate power law exponent 
r was measured for P(s, /) with s < 2 14 ~ 1.6 x 10 4 . To 
improve the precision we can use a very accurate value for 
T r = 1 + (t— l)//x = 1.299(3) measured by direct numeri- 
cal integration of Eq. (§) . Our resources allowed us to in- 
tegrate this equation forward for r < 2 14 , which is equiv- 
alent to measuring P(s,f) up to s = r 1 ^ ~ 1.5 x 10 7 , 
i.e. over much wider range than from numerical integra- 
tion of Eq. ([j]). Using exponent relations r = 2/i and 
t = 1 + /i(r r — 1) from t t = 1.299(3) one gets 



r = 1.176(2) 
H = 0.588(1) 



(9) 
(10) 



These are our best estimates of two basic exponents of 
the one-dimensional anisotropic BS model. The Monte 
Carlo simulations of the anisotropic BS model in d = 1 
are in perfect agreement with these values for r and /x. 
Indeed, Head and Rodgers |TJ| found pb = 0.59(3), while 
in Q it was measured to be fj, — 0.60(1). In Q they 
also measured the exponent tt — 2.42(5) of the distribu- 
tion of spatial jumps of the minimal site. This should 
be compared to our prediction for this exponent based 
on the scaling relation n = 1 + (2 — r)//x = 2.401(6). 
We also have performed Monte Carlo simulations of the 
anisotropic BS model in one and two dimensions. In ID 



we found = 0.58(1) and t\d — 1.17(1) in agreement 
with [jl2 13 j, our analytical results, and the direct simula- 
tion of Eq. (0) . In two dimensions our Monte Carlo simu- 
lations give H2D = 0.83(1) and T2D = 1-35(1). As shown 
in Fig. 2 these exponents lie on the t(/x) curve, valid 
for a general BS model. However, in the 2D anisotropic 
BS model we don't know the exact value of 7 to fix the 
position of the exponents on this curve. 

We also tested our theoretical prediction (r) ~ 2/ (f c — 
f) (which led us to 7 = 2) against the direct numerical 
integration of Eq. (0) . We measured the first moment (r) 
of the distribution Q(r, /), obtained as described above 
by direct numerical integration of Eq. (|§|) . As expected 
the power law divergence of (r) both above and below 
f c has the exponent — 1. It can be clearly seen when 
numerical derivative d((r)~ 1 )/df is plotted as a function 
of / (see Fig. 3). This derivative approaches different 
finite limits as / — > f c ± 0. The divergence of (r) at 
f c is clearly cut off by the finite size R of avalanches 
considered. This is illustrated in Fig. 3 by showing two 
curves for two different cutoff sizes B = 2 12 and 2 . 
From the asymptotical value of d((r)~ 1 )/df as / — > f c — 
we measure 1/7 = 0.498(3), which yields 7 = 2.01(2), in 
complete agreement with Eq. (§h. 

We were able to estimate yet another critical exponent 
from this plot. As it was shown in fl, the divergence of 
(r) in the overcritical regime above f c is given by (r) = 
(3/(f — f c ), where (3 = (r — 1)/(1 + fi — r) is the "order pa- 
rameter" exponent describing the scaling of the probabil- 
ity to start an infinite avalanche Poo(f) ~ (/ — jc)' 3 - For 
Poo = 1 - Y,7=i p ( s ' /) tnc Ec L- (§) readily gives dfPoo = 
(s^lPoo- Therefore, above the critical point one has 
(r) = {si*) - 0/Onpoo) = f3/(f~fc) + 0((f-f c r 2 )l The 
quadratic fit to d((r)~ 1 )/df above f c gives 1/(3 = 2.34(2), 
or (3 — 0.427(4). This is in excellent agreement with our 
theoretical prediction (3 = (r - 1)/(1 + /i - r) = 0.427(6) 
based on the best estimate of r r = 1.299(3) and the ex- 
ponent relation (3 = (r r — l)/(2 — r r ). 

In summary, we have analyzed the behavior of the 
anisotropic Bak-Sneppen model. We demonstrated that 
a nontrivial relation between critical exponents r and 
fj,, recently derived for the isotropic Bak-Sneppen model 
H||, holds for its anisotropic version as well. The ex- 
ponents measured by Monte Carlo simulations of the 
anisotropic Bak-Sneppen model in one and two dimen- 
sions are in agreement with this relation. For one- 
dimensional anisotropic Bak-Sneppen model we derive 
a novel exact equation (|jj) for the distribution Q(r, f) 
of avalanche spatial sizes. We also propose analogous 
equations for one-dimensional isotropic BS model and 
anisotropic BS model in d > 1. By studying the be- 
havior of the first moment of the distribution Q(r, /) we 
managed to extract the exact value 7 = 2 for one of the 
critical exponents of the one-dimensional anisotropic BS 
model. The values of critical exponents r and fi = d/D 
were found as coordinates of the intersection point be- 
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tween r(/i) and j(r,fx) = 2 curves. They are in excel- 
lent agreement with both Monte Carlo simulations of the 
model as well as results of numerical integration of the 
master equation for Q(r, /). We summarize our best es- 
timates for the exponents in one and two dimensions in 
Table 1. 

One of us (S.M.) would like to thank N.D.Mermin for 
asking the question that triggered this study. The work 
at Brookhaven National Laboratory was supported by 
the U.S. Department of Energy Division of Material Sci- 
ence, under contract DE-AC02-98CH10886. 



Exponent ID anisotropic BS 2D anisotropic BS 

t 1.176(2) 1.35(1) 

(j, 0.588(1) 0.83(1) 

a 0.412(1) 0.48(2) 

7 2 1.35(4) 

/3 0.427(6) 0.73(4) 

D 1.701(3) 2.41(3) 

7T 2.401(6) 2.57(3) 

TABLE I. The results from this table are obtained by nu- 
merical integration of Eq. ([]) for d — 1 and Monte Carlo 
simulations for d = 2. The existing exact exponent relations 
were then applied. 
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FIG. 1. The effective power law exponent >(R), measured as -r,r ff, (i?) = (log Q(R, /)-log Q(R-1, /))/(log_R-log(.R-l)) 
for two different /. Since our method is free from finite size effects, one can be sure that f c is in (1.286, 1.287) and r r = 1.299(3). 
Q(R, f) was obtained by numerical integration of Q with R < -R m ax = 2 14 = 16384. A second order Runge Kutta method 
with 5f = 10 -3 was used. 
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FIG. 2. The curve r(/n) from Ref. 9 is shown along with the results of Monte Carlo simulations of both isotropic (o) and 
anisotropic (•) BS models in one and two dimensions. The intersection point of r(/i) with the line r = 2/i determines the 
exponents of the ID anisotropic BS model in agreement with Monte Carlo results. 



7 




8 



